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EXECUTIVE  SUMMARY 


Large-scale  molecular-dynamics  (MD)  simulations  were  performed  to  investigate:  i)  sintering 
process,  structural  correlations,  and  mechanical  behavior  including  dynamic  fracture  in 
microporous  and  nanophase  Si3N4;  ii)  crack-front  propagation  in  amorphous  silica;  iii) 
hypervelocity  impact  damage  in  diamond  films  and  crack-front  instabilities  in  graphite;  iv)  phonons 
in  graphitic  tubules;  and  v)  amorphization  and  fracture  in  nanowires.  The  simulations  were  carried 
out  with  highly  efficient  multiscale  algorithms  and  dynamic  load-balancing  schemes  for  mapping 
irregular  atomistic  simulations  on  distributed-memory  parallel  architectures.  These  research 
activities  resulted  in  twenty-four  journal  publications  and  thirty  invited  presentations.  Two 
graduate  students  received  dual-degrees— Ph.D.  in  physics  and  a  M.S.  from  the  Department  of 
Computer  Science— during  the  course  of  this  project. 

This  project  involved  two  faculty  members,  two  graduate  students,  and  a  postdoc.  It  also  involved 
interactions/transitions  with  researchers  at  governmental  laboratories  and  industry. 
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BACKGROUND 


Advanced  ceramics  are  excellent  materials  for  applications  requiring  extreme  operating  conditions. 
The  capability  to  withstand  high  temperatures,  combined  with  their  high  strength  and  low  weight, 
makes  them  highly  desirable  for  aerospace  applications  such  as  high  thrust-to-weight  ratio  turbine 
engines,  hypersonic  aerospace  vehicles,  and  space  satellites.  These  materials  are  also  very  useful 
for  surface  transportation,  electronics,  and  advanced  manufacturing  industries. 

The  fundamental  issue  concerning  ceramics  is  the  relationship  between  the  microscopic 
structure,  atomic  diffusion,  and  mechanical  and  thermal  properties.  An  understanding  of  this 
relationship  on  the  atomic  scale  is  critical  for  the  design  of  novel  ceramics  with  application-specific 
properties. 

This  DEPSCoR  project  involved  the  investigation  of  properties  and  processes  in 
nanostructured  ceramics.  Using  large  scale  molecular-dynamics  (MD)  simulations,  we 
investigated: 

•  Structural  correlations  and  mechanical  behavior,  including  crack-front  propagation,  in 
crystalline  and  nanophase  silicon  nitride; 

•  Growth  of  pore  interfaces  and  brittle  and  ductile  behavior  in  amorphous  silica; 

•  Mechanical  properties,  thermal  transport,  and  low-frequency  floppy  modes  in  microporous 
silicon  nitride; 

•  Early  stages  of  sintering  of  ceramic  nanoclusters; 

•  Damage  in  diamond  coatings  due  to  particle  impact; 

•  Crack  propagation  in  graphite; 

•  Strain  effects  and  phonons  in  graphitic  tubules; 

•  Amorphization  and  fracture  in  nanowires. 

These  MD  simulations  were  executed  with  highly  efficient,  portable  and  scalable, 
multiresolution  algorithms  including  the  fast  multipole  method  for  the  long-range  Coulomb 
interaction  and  a  dynamic  load-balancing  scheme  for  mapping  irregular  applications  on  parallel 
machines.  These  research  activities  have  resulted  in  24  journal  papers  and  30  invited  presentations. 
The  next  section  provides  an  overview  of  our  research  accomplishments. 

§2  RESEARCH  ACCOMPLISHMENTS 

§2.1  Structure  and  Mechanical  Properties  of  Nanophase  Si3N4 

Advanced  ceramics  such  as  Si3N4,  SiC,  Si02,  and  AljOj  are  highly  desirable  materials  for 
applications  requiring  extreme  operating  conditions.  Light  weight,  elevated  melting  temperatures, 
high  strengths,  and  wear  and  corrosion  resistance  make  them  very  attractive  for  applications  in 
aeronautics,  automotive,  electronics,  and  advanced  manufacturing  industries.  The  only  serious 
drawback  of  ceramics  is  that  they  are  brittle  at  low  to  moderately  high  temperatures. 

In  recent  years  a  great  deal  of  progress  has  been  made  in  the  synthesis  of  ceramics  that  are 
much  more  ductile  than  conventional  coarse-grained  materials.  (These  so  called  nanophase 
materials  are  fabricated  by  in-situ  consolidation  of  nanometer  size  clusters.)  Despite  a  great  deal  of 
research,  many  perplexing  questions  concerning  nanophase  ceramics  remain  unanswered. 
Experiments  have  yet  to  provide  information  regarding  the  morphology  of  pores  or  the  structure 
and  dynamics  of  atoms  in  nanophase  ceramics.  As  far  as  modeling  is  concerned,  only  a  few 
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aiomisuc  simulations  of  nanophase  materials  have  been  reponed  thus  far.  (This  is  due  to  the  fact 
that  these  simulations  are  computationally  very  demanding:  A  realistic  simulation  of  a  nanophase 
solid  requires  l()'’-l(f  time  steps  for  proce.ssing  and  ~  10^  atoms,  since  each  nanocluster  itself 
consisLS  of  KK-lO'*  atoms.) 

U.sing  1.08  million  particle  MD  simulations  we  have  investigated  the  structure  and  mechanical 
behavior  of  nanophasc  Si3N4.  The  MD  simulations  are  performed  with  effective  two-  and  three- 
body  interatomic  potentials:  The  two-body  terms  consist  of  steric  repulsion  between  atoms, 
screened  Coulomb  potentials  due  to  charge  transfer  between  Si  and  N,  and  a  charge-dipole 
interaction  that  takes  into  account  the  large  electronic  polarizability  of  nitrogen;  the  three-body 
potential  in  Si3N4  takes  into  account  covalent  effects  through  bond-bending  and  bond-stretching 
terms.  To  establish  the  validity  of  this  interaction  scheme,  we  compare  the  MD  results  with  a 
variety  of  experimental  measurements.  We  fmd:  i)  the  bond  lengths  and  bond-angle  distributions 
in  crystalline  and  amorphous  Si3N4  are  in  excellent  agreement  with  experiments;  ii)  the  positions 
and  relative  heights  of  the  peaks  in  the  static  structure  factor  for  amorphous  Si3N4  are  in  good 
agreement  with  neutron  scattering  measurements:  iii)  the  phonon  density-of-states  of  crystalline  a- 
Si3N4  agrees  well  with  inelastic  neutron  scattering  experiments;  iv)  the  specific  heat  of  crystalline 
a-Si3N4  is  in  excellent  agreement  with  experiments  over  a  wide  range  of  temperatures;  and  v)  the 
bulk  modulus  and  Young's  moduli  along  different  directions  in  a-SijN^  deviate  less  than  10% 
from  the  experimental  values. 


Figure  1:  Snapshot  of  a 
configuration  of  nanophasc 
SiiNj  before  consi)lidatiiHi 
(left).  The  other  picture  shows 
the  system  consolidated  witli 
an  applied  pressure  of  15  GPa. 
The  densitN  of  the 
consolidated  sNsieiu  is  2,94 
g/cc  which  is  92''(  of  the 
density  of  the  (x-Si,N,  crystal. 


Nanophase  Si3N4  is  generated  by  consolidating  a  random  cluster  configuration  (108  clusters 
with  10,052  atoms  each)  with  the  constant-pressure  MD  approach.  Figure  1  shows  the  nanophase 
systems  before  and  after  consolidation  under  an  applied  pressure  of  15  GPa.  Pair-distribution 
functions  and  bond-angle  distributions  reveal  that  inlerfacial  regions  in  the  consolidated  nanophase 
Si3N4  arc  amorphous  and  they  contain  a  large  number  of  undercoordinated  Si  atoms.  Systems 
sintered  at  low  pressures  (~  1  GPa)  have  percolating  pores  whose  surface  roughness  exponents 
(0.46  and  0.86)  are  in  excellent  agreement  with  experiments.  We  also  find  that  the  dependence  of 
elastic  moduli  on  porosity  and  grain  size  in  nanophase  Si3N4can  be  understood  in  terms  of  a  three- 
phase  model  for  heterogeneous  materials  (please  see  Fig.  2). 
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Porosity 


Fisjure  2:  Porosil\  dopendeiKe  ul  ihc 
hulk  modulus  iKi  and  shear  modulus  (G), 
Open  circles  denote  Ihe  MD  resulls  lor  ihe 
hulk  modulus  of  nanophase  silicon  nitride 
with  cluster  si/e  D=h()A:  open  iriantiles 
and  squares  represent  the  MD  resulls  lor 
hulk  and  sheiu'  moduli,  respectively,  of  the 
nanophase  s\siem  with  45A  clusters. 
Solid  circles  and  iriuni^les  denote  hulk 
moduli  and  svdid  squares  represent  shear 
moduli  calculated  livtm  the 
multicomponent  model  |Budiansk\.  J, 
Mech,  Plus,  Solids  13.  223  (1965)1.'  The 
solid,  dashed,  and  dotted  lines  are  drawn  to 
guide  the  eye. 


§2.2  Cr.4CK-Front  PROP.4G.\rioN  IN  Nanoph.\se  Silicon  Nitride 

Undersianding  the  role  of  microstruciures  in  fraciure  i$  one  of  the  most  challenging  problems  in 
materials  science.  Nanophase  materials  are  ideal  systems  to  examine  this  issue  at  the  atomistic 
level  since  microstructures  in  these  materials  are  only  a  few  nanometers  in  dimension.  Using  1.08 
million  particle  MD  simulations,  we  have  investigated  the  morphological  and  dynamic  aspects  of 
fracture  in  nanophase  Si3N4.  Figure  3  shows  snapshots  of  the  crack  front  at  various  values  of  the 
external  strain.  The  crack  front  consists  of  pores  connected  to  the  notch  (shown  in  magenta;  the 
remaining  disconnected  pores  are  shown  in  orange).  Figure  3  (a)  shows  a  snapshot  taken  10  ps 
after  the  notch  is  insened.  We  observe  initial  development  of  the  crack  front  and  the  growth  of  a 
few  crack  branches  in  the  system.  These  local  branches  and  nanoclusters  tend  to  arrest  the  motion 
of  the  crack  front,  and  further  crack  propagation  is  only  possible  if  additional  strain  is  applied. 
Figure  3(b)  shows  a  snapshot  of  the  crack  from  in  the  nanophase  system  under  11%  strain. 
Comparing  with  Fig.  3(a),  it  is  evident  that  the  crack  from  has  advanced  significantly  and 
coalesced  with  the  pores  in  the  center.  We  also  observe  tlrat  pores  and  inierfacial  regions  allow  die 
crack  front  to  meiinder  and  form  a  branched  structure.  Figure  3(c)  shows  a  crack-front  snapshot 
10  ps  after  the  system  was  strained  to  14%'.  At  this  time  a  secondary  crack  (top  left  hand  comer  of 
the  figure)  with  several  local  branches  merges  widi  the  primary  crack.  With  further  increase  in  the 
strain  the  secondary  from  advances  toward  the  initial  notch  while  the  crack  keeps  growing  laterally. 
When  the  strain  reaches  30%>,  the  crack  finally  separates  the  material  into  two  disconnected  parts 
and  the  system  is  now  completely  fractured.  It  should  be  noted  that  the  critical  strain  (30%)  at 
which  the  nanophase  system  fractures  is  enomious  compared  to  what  the  crystalline  Si3N4  system 
can  sustain  (at  an  applied  strain  of  only  3%  the  crystal  undergoes  cleavage  fraemre).  This  is  due 
to;  i)  plastic  defomiation  in  interfacial  regions;  ii)  deflection  and  arrest  of  cracks  by  nanoclusters; 
and  iii)  multiple  crack  branching.  None  of  these  mechanisms  are  operative  in  the  crystalline 
system.  This  demonstrates  the  dramatic  effect  of  nanostructures  on  fracture. 

To  compare  the  toughness  of  nanophase  and  crystalline  systems,  we  have  calculated  the 
fracture  energy  per  unit  area  in  both  nanophase  and  crystalline  Si3N4  systems  (with  the  same 
geometry).  For  the  ( 10  T  0)  surl'ace  of  an  atomically  sharp  mode  I  crack  (film  surlhce  was  ((K)()  1 )) 
in  an  a-Si3N4  crystalline  film,  the  stress  distribution  was  calculated  using  the  MD  method.  From 
the  prefacior  of  the  tcmi  of  the  local  stress  (r  is  the  distance  from  the  crack  tip),  we  have 
estimated  the  critical  surcss-imensity  factor.  This  analysis  gives  a  value  which  is  close  to  our  MD 
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estimate  of  tlte  intrinsic  toughness  (1.3  MPa*m'’  ^).  From  the  calculations  of  the  fracture  energy  per 
unit  area,  we  estimate  that  the  toughness  of  the  nanophase  Si3N4  is  six  times  larger  than  that  of  the 
crystal. 


(a) 

_  (bf 

_ ,  __  (cL 

iy 

; 

1 

■f 
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Figure  3:  Snapshots  of  the  crack  front  (shown  in  magenta)  along  witli  large  (>6.4  nni')  isolated 
pc)res  (shown  in  orange)  in  the  strained  nanophase  system:  a)  liritial  notch  with  crack  branches  and 
pores  in  the  system  under  an  applied  strain  of  b)  the  primary  crack  front  after  the  strain  is 
increased  to  \\’/r.  and  c)  the  primary  and  secondary  crack  fronts  at  14'/f  strain  on  the  system. _ 


We  have  also  investigated  the  morphology  of  crack  fronts  in  the  fractured  nanophase  Si3N4.  In 
recent  years  this  issue  has  drawn  a  great  of  experimental  attention.  Measurements  on  a  variety  of 
brittle  and  ductile  solids  have  revealed  that  the  widths  of  fracture  profiles  scale  with  the  lengths  as 
H'  ~  Lr.  For  out-of-plane  fracture  profiles,  the  roughness  exponent,  above  a  certain  length 
scale,  has  a  value  of  0.8  in  three  and  0.7  in  two  spatial  dimensions.  Many  experiments  indicate 
that  these  values  of  (j.  may  be  independent  of  the  material  or  the  mode  of  fracture.  At  small  length 
scales  or  when  the  crack  front  propagates  quasi  statically,  the  roughness  exponent  crosses  over  to  a 
smaller  value  ~  0.5.  Recendy  it  has  also  been  pointed  out  that  there  is  another  distinct  roughness 
exponent,  ^  =  0.6,  associated  with  in-plane  fracture  profiles. 


Figure  4:  i.d  lug-log  ploi  of  ihc  iK'ighl-bcighi  coriclalion  limclion.  versus  :  for  ihe  oul-ot- 

plane  liLiciure  profile  o.-j;  (b)  ihe  van.iiion  ol  lor  (he  in-plaiie  Iraclure  profile  \i:i. 


To  investigate  the  nature  of  self-affme  fracture  surfaces  in  nanophase  Si3N4,  we  calculate  the 
height-height  correlation  functions  both  in  and  out  of  the  fracture  plane  y-c.  Figure  4a  shows  that 
the  best  fit  to  the  out-of-plane  height-height  correlation  function  g,^(z)  (=  <[x(z+z^j)  -  xizo)]'>‘^) 

for  the  fracture  profile  x(z)  requires  two  roughness  exponents:  ^  =  0.84±0. 12  above  a  certain 
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length  scale  (64A)  and  -  ().58±().  14  otherwise.  For  the  other  out-of-plane  height-height 
correlation  function  gxx(y)^  the  roughness  exponent  Q|  =  0.75±0.08.  The  MD  results  for  ^j.  and  Q| 
are  very  close  to  experimental  values.  We  have  also  determined  the  roughness  exponent,  for  the 
in-plane  fracture  profile  y(z):  Figure  4b  shows  that  the  best  fit  to  the  corresponding  height-height 
correlation  function  gives  ^  =  ().57±().08.  Measurements  of  in  Al-Li  and  Superai  Ti3Al-based 
alloys  yield  0.60±0.04  and  0.54±().()3,  respectively  [E.  Bouchaud,  J.  Phys.  Cond.  Matter]. 

§2.3  Sintering  of  Silicon  Nitride  N.vnoclusters 


The  study  of  structural  and 
mechanical  properties  of  nanophase 
SijN^  was  preceded  by  MD 
simulations  of  sintering  of  silicon 
nitride  nanoclusters  (each  cluster 
consisting  of  20,335  atoms).  Our 
MD  calculations  provide  a 
microscopic  view  of  neck  formation 
during  early  stages  of  sintering 
(Fig.  5).  In  the  case  of  Si3N4 
nanocrystals  at  2000K.  we  initially 
observe  considerable  relative 
motion  of  clusters.  Subsequently  a 
few  Si  and  N  atoms  join  the  two 
nanocrystals  and,  thus  bound,  they 
continue  to  rotate  relative  to  each 
other  for  100  ps.  In  the  next  1(X)  ps,  the  relative  motion  subsides  and  we  observe  a  steady  growth 
of  an  asymmetric  neck  between  the  two  nanocrystals.  In  the  neck  region,  there  are  more  four-fold 
than  three-fold  coordinated  Si  atoms.  We  have  also  investigated  the  sintering  of  amorphous  sUicon 
ninide  nanoclusters  at  2,000K.  The  neck  between  amorphous  nanoclusters  is  much  more 
symmetric  tlian  the  neck  between  themiaUy  rough  nanocrystals.  The  neck  region  between 
amorphous  nanoclusters  has  nearly  the  same  number  of  three-  and  four-fold  coordinated  Si  atoms. 
For  both  nanocrystals  and  amorphous  nanoclusters,  sintering  is  driven  by  diffusion  of  surface 
atoms.  The  diffusion  in  the  neck  region  of  amorphous  clusters  is  four  times  faster  than  in  the  neck 
between  nanocrystals. 

§2.4  Floppy  Modes  and  Therivial  Transport  in  Microporous  Si3N4 

Physical  properties  of  Si3N4  are  greatly  influenced  by  voids.  (They  are  usually  present  in  the 
system  as  a  result  of  synthesis  by  chemical  vapor  deposition  or  other  processes.)  Using  MD 
simulations  we  have  investigated  the  elTect  of  porosity  on  mechanical  properties  and  themial 
conductivity  of  amorphous  Si3N4.  Under  uniform  dilation,  pores  begin  to  fomt  at  a  mass  density 
of  p  =  2.6  g/cc.  With  further  decrease  in  the  density,  the  number  of  pores  and  their  sizes  increase 
dramatically  and  around  2.0  g/cc  the  largest  pore  percolates  through  the  system.  Near  percolation 
the  average  pore  volume  diverges  as  (p  -  p^,)”'  with  y  =1.95±U.  17.  This  is  consistent  with  the 
percolation  theory. 

The  themial  conductivity  of  microporous  Si3-N4  is  calculated  from  the  heat  current-current 
correlation  function  and  a  non-equilibrium  MD  approach.  Between  300  and  1500K,  the  themial 
conductivity  scales  as  p^  where  r  =1.5.  Experimental  data  on  carbon  and  silica  aerogels  are  in 


Figure  5:  Snapshots  showing  sintering  of  Si,N,  nanodusters 
at  2,()()()K:  (a)  nanocrystals  at  time  l  =  O;  (b)  sintered 
nanocrysials  at  t  =  200  ps;  and  (c)  sintered  amorplious 
nanodusiers  at  t  =  400  ps. 
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agreement  with  this  result.  The  MD  results  for  Young’s  modulus  E  vary  as  with  T  =3.6. 
This  is  also  close  to  the  experimental  value  of  T  for  silica  aerogels. 

The  nature  of  low-energy  floppy  modes  in  microporous  Si3N4  has  also  been  investigated  with 
the  MD  approach.  Floppy  modes  appear  continuously  in  amorphous  Si3N4  as  the  network 
connectivity  is  reduced;  in  ^e  crystal  they  appear  suddenly  at  20%  volume  expansion.  We  find 
that  their  contribution  to  density-of-states  varies  linearly  with  energy  and  they  significantly  enhance 
the  specific  heat  of  microporous  Si3N4. 

§2.5  Hyper  VELOCITY  Impact  Damage  in  a  Diamond  Coating 

Recently  researchers  at  Phillips  Lab’s  Space  Environment  Interaction  Branch  have  designed  an 
experimental  technique  to  simulate  debris  impact  in  low  earth  orbit  (LEO)  using  micro-flyer  plates. 
The  latter  simulate  debris  particles  traveling  at  speeds  in  the  range  of  8-10  km/s.  The  plates  are 
propelled  by  a  laser  beam  inside  a  space  chamber.  Currently  the  debris  experiments  are  performed 
with  an  expensive  60-foot-long  gun. 

In  an  attempt  to  simulate  such  debris-impact  experiments  on  parallel  computers,  we  have 
performed  large-scale  MD  calculations  to  determine  the  hypervelocity  impact  damage  caused  by  a 
diamond  crystallite  on  a  diamond  coating.  (Rocketdyne  in  collaboration  with  Phillips  Laboratory 
are  using  chemical  vapor  deposition  to  apply  thin  films  of  diamond  to  surfaces  for  a  wide  range  of 
rocket  applications.)  The  MD  simulations  are  based  on  a  reactive  bond  order  (REBO)  potential 
proposed  by  Brenner  [Phys.  Rev.  B  42,  9458  (1990)].  This  potential  is  quite  unique  in  that  it 
treats  covalent  bonding  in  molecular  and  solid-state  structures  with  a  single  classical  expression. 
Moreover,  the  potential  can  describe  chemical  processes  such  as  covalent  bond  formation  and 
breaking.  The  parameters  of  the  potential  have  been  fitted  to  experimental  values  and  first- 
principles  calculations  of  bond  lengths,  bond  energies,  and  force  constants  for  several  solid-state 
and  molecular  systems.  The  REBO  potential  has  been  shown  to  provide  a  good  description  of 
various  properties  outside  the  fitting  database  such  as  elastic  constants  and  vibrational  spectra  of 
diamond  and  in-plane  elastic  properties  of  graphite.  The  results  for  the  surface  reconstruction  and 
energetics  of  point  defects  in  diamond  are  also  in  good  agreement  with  experiments  and  electronic- 
structure  calculations.  The  calculated  formation  energies  for  an  extensive  set  of  hydrocarbon 
molecules  agree  well  with  experimental  results. 

We  have  performed  three  hypervelocity  impact-damage  simulations.  In  the  first  case  the 
crystallite  strikes  the  diamond  film  with  a  velocity  of  8  km/s.  The  crystallite  penetrates  the  film  to  a 
depth  of  20  A,  causing  the  diamond  structure  of  the  film  to  graphitize  without  any  crater  formation. 
Subsequently  the  nearly-intact  crystallite  emerges  from  the  coating  at  a  speed  of  4  km/s.  In  the 
second  simulation  the  incoming  crystallite  strikes  the  coating  at  a  speed  of  1 1  km/s.  Part  of  the 
crystalhte  disintegrates  after  impact,  and  almost  all  of  its  kinetic  energy  is  imparted  to  the  coating. 
In  this  case  we  observe  the  formation  of  a  crater  along  with  amorphization  in  the  impact  region.  In 
the  third  case  the  crystallite  hits  the  diamond  coating  at  a  velocity  of  15  km/s.  Unlike  the  first  two 
cases,  the  crystallite  completely  disintegrates  causing  extensive  damage  including  the  creation  of  a 
crater. 
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Filjure  6:  The  first  nine  snapshots  show  the  impact  of  a  cliamond  crystallite  on  a  diamond  coating. 
I'hese  are  the  results  of  our  million  atom  MD  simulations  over  a  time  period  of  28  pico  .seconds. 
The  crystallite  hits  the  coating  at  a  speed  of  8  km/s.  Alter  impact,  it  penetrates  the  coating  and  then 
hounces  out  at  a  speed  of  4  km/s,  Tlie  last  snapshot  shows  that  the  energy  imparted  to  the  coaling 
causes  graphiti/ation  under  the  coniaci  region. 


§2.6  Dynamic  Fracture  in  a  Graphite  Sheet 

Recently  we  have  also  investigated  crack  propagation  in  a  graphite  sheet  using  million-atom  MD 
simulations  on  parallel  machines.  (These  simulations  are  also  based  on  the  reactive  bond  order 
potential  of  Brenner.)  Fracture  dynamics  is  investigated  under  a  constant  applied  strain  after 
inserting  a  notch  of  length  20A.  We  consider  two  different  orientations  of  the  graphite  sheet, 
which  exhibit  completely  different  behavior  during  crack-front  propagation: 

•  Under  an  applied  uniaxial  strain  of  12%  the  system  G(l,l),  in  which  a  fraction  of 
covalent  bonds  are  parallel  to  the  strain  (Fig.  7b),  undergoes  cleavage-like  fracture; 

•  For  the  G(l,())  orientation,  the  crack  front  develops  multiple  secondary  branches  and 
overhangs  (Fig.  7a)  when  its  speed  reaches  approximately  half  the  Rayleigh  wave 
speed.  Within  the  same  secondary  branch  the  crack-front  profile  has  a  roughness 
exponent  ^  =  0.5,  whereas  for  interbranch  fracture  surface  profiles  ^  =  0.7.  The 
scenario  of  multiple  kx'a!  branches  and  the  values  of  the  roughness  exponenLs  are  in 
agreement  with  recent  exjjcrimental  investigations  of  fracture. 
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We  have  also  determined  the 
fracture  toughness  of  the 
graphite  sheet  in  the  G(l.l) 
orientation  using  the  MD 
approach.  The  system  with  a 
notch  of  50  A  was  subjected 
to  different  strains. 
Applying  the  Griffith 
criterion,  with  the  surface 
energy  determined  from  the 
final  fracture  surface,  the 
critical  stress  is  estimated  to 
be  7 1  GPa.  The  crack  starts 
propagating  at  a  stress  of  66 
GPa,  which  corresponds  to  a 
critical  strain  of  10%.  From 
these  results,  the  value  of  the 
fracture  toughness  (defined 


as  (T/rVc ,  where  cr^  is  the  critical  stress  at  which  the  crack  starts  propagating  and  c  is  the  length 


of  the  notch)  is  found  to  be  4.7  MPa»m'^^  We  have  also  calculated  the  local-stress  distribution 


around  the  notch  just  before  the  onset  of  crack  propagation.  The  data  are  well-described  by  1/ r 


dependence  with  a  stress-intensity  factor  of  6  MPa»m‘'l  The  angular  dependence  of  the  local 


stresses  is  also  in  good  agreement  with  the  elasticity  theory. 


§2.7  Amorphization  and  Fracture  in  Nanowires 

Silicon  diselenide  is  an  important  chalcogenide  material.  It  is  one  of  the  basic  solid  electrolytes  for 
advanced  high-energy  electrochemical  cells.  Recently  we  investigated  the  dynamics  of 
amorphization  and  fracture  in  SiScj  nanowires  using  large-scale  MD  simulations.  The  calculations 
are  based  on  interatomic  potendals  which  combine  charge-transfer,  electronic  polarizabilities,  steric 
effects,  and  three-body  covalent  interactions.  The  validity  of  the  interaction  scheme  is  established 
by  comparing  the  MD  results  with  various  available  measurements.  We  find:  structural 
correlations  in  both  crystalline  and  amorphous  systems  are  in  excellent  agreement  with 
experiments;  the  calculated  melting  temperature  (1250+2UK)  is  very  close  to  the  experimental  value 
(1233±5K);  and  the  vibrational  spectra  of  crystalline  and  amorphous  SiSCj  are  also  in  good 
agreement  with  neutron  scattering  measurements. 


The  MD  simulations  for  SiSCj  nanowires  reveal  that  the  cTitical  strain  (~  15%)  of  fracture  i.s 
independent  of  the  number  of  chains  in  the  nanowire.  Also  fracture  is  preceded  by  amorphization 
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in  nanowires.  At  first,  bonds  between  edge-shaiing  tetrahedra  break  in  the  chains  at  the  outennosi 
layer.  These  broken  chains  form  comer-sharing  tetrahedra  by  cross-linking  with  neighboring 
chains.  Lxx'al  amorphization  propagates  across  the  nanowire  to  form  a  sandwich-like  structure 
with  an  amorphous  region  in  the  middle.  Cracks  then  nucleate  at  the  boundaries  and  propagate 
inward  from  the  nanowire  surface  until  the  nanowire  fractures  (Fig.  8  ). 

§2.8  Structural  Correlatio.ns  and  BRiTTUE/DucrtLE  Behavior  in  Amorphous 

Silica 

Silica  has  numerous  technological  applications  in  microelectronic  devices  and  micro-electro¬ 
mechanical  systems  (MEMS).  Atomistic  simulations  have  been  used  to  study  the  behavior  of 
silica.  We  have  investigated  the  fracture  process  in  vitreous  SiO,  by  applying  the  parallel 
multiresolution  MD  approach  to  a  system  of  1.18  milhon  atoms.  The  calculations  are  based  on  an 
interatomic  potential  which  includes  the  effects  of  charge  transfer,  electronic  polarizabilities,  and 
steric  repulsion  through  two-body  temis;  covalent  effects  are  included  via  three-body  bond¬ 
bending  and  bond-stretching  temis.  Figure  9  shows  a  comparison  between  the  MD  and  neutron- 
scattering  results  for  the  static  structure  factor,  S(q),  of  amorphous  silica.  Clearly  the  simulation 
results  for  S(q)  are  in  very  good  agreement  with  experiments.  The  calculations  for  elastic  nioduh 
and  the  phonon  density-of-states  of  amorphous  SiO,  are  also  in  good  agreement  with  experimental 
measurements. 


Figure  9:  Neuiron-scauering  striicluro 
I'aclor  of  aniorplious  SiO>.  The  MD 
resiiUs  agree  well  with  those  Ifuin 
neutron  diffraction  experiments 
|.lohnson,  Wrielil.  and  .Sinclair.  .1.  Non- 
Cryst.  Solids  58.  100  (iy.S.M|. 

including  the  height  ol  the  first  sharp 
diffraction  peak 
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Figure  10  shows  the  MD  results  for  fracture  (snapshots  of  crack  fronts)  in  amorphous  sihea  at 
300K  and  900K.  It  is  evident  that  the  system  undergoes  cleavage  fracture  at  300K.  At  the  critical 
strain  for  fracture  (8%),  the  crack  tip  reaches  a  terminal  speed  of  1.40  km/s  which  is  nearly  half  the 
Rayleigh  wave  speed.  (At  600K,  the  fracture  is  again  cleavage-hke  with  the  temiinal  crack-tip 
speed  reaching  a  value  of  0.85  km/s.)  In  contrast,  at  9(X)K  the  amorphous  silica  system  shows 
plastic  behavior  up  to  a  strain  of  17%.  Currently  we  are  investigating  sintering,  suaictural 
correlations,  and  mechanical  properties  of  nanophase  amorphous  silica  .system. 
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Figure  10:  Snapshots  showing  brittle  fracture  in  vitreous  silica  at  3(X)K  (the  top  three  figures).  The 
boitoni  three  figures  .show  ductile  behavior  at  9()()K,  Only  regions  of  depleted  density  are  shown. 

§2.9  Atomistic  Processes  in  Silicon  Carbide  and  Alumina 

Silicon  carbide  and  alumina  are  excellent  simcwral  ceramics  for  applications  requiring  extreme 
operating  conditions.  We  are  investigating  properties  and  processes  in  these  materials  using  large- 
scale  MD  simulations.  The  SiC  simulations  are  based  on  an  empirical  bond-order  potential 
propo.sed  by  Tersoff  [Phys.  Rev.  B  39,  5566  (1989)].  The  MD  resulLs  for  structural  properties 
(bond  lengths,  bond  angles,  etc.),  elastic  moduli,  and  the  cohesive  energy  of  crystalline  SiC  are  in 
exceUem  agreement  with  experiments.  The  melting  temperature  obtained  from  MD  simulations 
using  Tersoff  s  potential  agrees  well  with  the  result  of  the  electronic-structure  calculations. 
Currently  we  are  perlhrming  MD  simulations  to  understand  the  crack-front  morphology  and  the 
effect  of  strain  rate  on  dynamic  fracture  in  crystalline  and  amorphous  SiC  systems  Preliminary 
results  of  these  fracture  simulations  are  shown  in  Fig.  11. 

The  MD  simulations  for  aluminum  oxide  are  based  on  an  interaction  scheme  proposed  by 
Streitz  and  Mintmire  [Phys.  Rev.  B  50.  11996  (1994)].  This  potential  allows  the  local  atomic 
charges  to  vary  with  the  environment.  The  dynamic  charge  transfer  and  the  resulting  electrostatic 
interaction  is  merged  with  an  embedded- atom  potential  yielding  a  scheme  which  can 
simultaneously  treat  FCC  aluminum  and  a-alumina  systems.  This  potential  has  been  optimized  to 
reproduce  the  bulk  structural  and  mechanical  properties  of  FCC  aluminum  and  a-alumina. 
Additionally,  the  potential  yields  reasonable  values  for  surface  energies  and  surface  relaxations. 

We  have  implemented  the  MD  simulations  for  Al^O^  using  an  extended  Lagrangian  scheme  in 
which  the  charges  are  treated  as  dynamical  variables.  The  long-range  Coulomb  interaction  is 
implemented  with  an  0{N)  fast  multipole  method  (FMM).  We  have  already  determined  the  phonon 
density-of-states  (DOS),  infrared  spectrum,  and  Raman  spectrum  of  a-A]203  crystal.  The  MD 
results  compare  well  with  experimental  measurements.  Currently,  we  are  investigating  the 
behavior  of  Al/AljO,  interface.  The  preliminary  results  are  shown  in  Fig.  1 1 . 
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Figure  11.  The  figure  on  the  left  shows  the  unset  ul  fracture  in  a  million  atom  MD  sitnulaiion  of 
crystalline  SiC  under  !'/(  strain  at  ?()()K.  The  other  figure  shows  dynamic  charge  transfer  at 
Al/AfO,  interface.  The  interface  invof  es  1 1  1  1 1  ,A1  and  ((tool )  .Al-O,  surfaces.  The  white  circles 
show  the  atomic  positions  (the  small  circle  is  aluminum  anil  the  large  circle  is  oxygen)  The 
color  bar  shows  the  valence  of  aluminum  atoms. 


§2.10  P.4RALLEL  Multiresolution  Molecular-Dynamics  Simulations 

Atomistic  simulations  discussed  in  §2.1  -  §2.9  require  considerable  computational  resources 
because  of  long  processing  times,  large  systems,  and  compute-intensive  interatomic  interactions 
(e.g.  long-range  Coulomb  and  three-body  covalent  forces).  These  simulations  are  not  feasible 
without  efficient  algorithms  or  massively  parallel  machines. 

We  have  designed  a  highly  efficient,  multiresolution  molecular-dynamics  (MRMD)  approach  to 
carry  out  large-scale  materials  simulations  on  parallel  machines.  This  approach  includes  the  fast 
multipole  method  (FMM)-a  divide-and-conquer  scheme  which  reduces  the  computational  cost  for 
the  Coulomb  interaction  from  0{hfi)  to  0(N).  The  FMM  is  based  on  a  hierarchical  subdivision  of 
space  (Fig.  12),  where  at  each  spatial  scale  the  influence  of  a  large  collection  of  charges  is 
approximated  to  any  desired  precision  by  a  multipole  expansion.  The  method  is  well-suited  to 
parallel  architectures;  Unlike  the  direct  method  which  must  communicate  the  locations  and 
strengths  of  a  large  number  of  charges  between  distant  processors,  the  FMM  transmits  a 
compressed  representation  of  this  information,  namely  the  coefficients  of  the  multipole  expansion. 
In  addition  to  FMM,  the  MRMD  approach  includes  a  multiple  time-scale  algorithm  for  shori- 
ranged,  two-body  interactions  and  a  separable  scheme  for  three-body  covalent  interactions. 


Figure  12:  Multilevels  in  the  FMM  for  a  2D  system  (left);  execution  (solid  curves)  and 
communication  times  (dashed  curves)  per  MD  time  step  for  the  MRMD  approach  (right),  p  is  the 
number  of  processors. 


Figure  12  shows  the  perfonnance  of  the  MRMD  approach  on  the  512-node  Touchstone  Delta 
machine  at  Caltech  and  the  128-processor  IBM  SP  system  at  Argonne  National  Laboratory.  The 
execution  time  scales  linearly  with  the  system  size  and  the  computation  time  dominates  the 
communication  time.  Using  MPI  (Message  Passing  Interface),  the  MRMD  has  been  ported  to 
other  parallel  platforms  as  well. 

§2.11  Dynamic  Load-Balancing  Scheme  for  Irregular  Applications 

In  many  atomistic  simulations  (e.g.  porous  materials  or  fracture)  one  encounters  a  non-uniform 
spatial  distribution  of  atoms.  Regular  partitioning  of  such  a  system  into  subsystems  of  equal 
volume  and  then  one-to-one  mapping  of  subsystems  on  processors  of  a  parallel  machine  leads  to 
an  imbalance  of  the  computational  workload  among  the  processors.  A  load-balancing  scheme  has 
been  designed  by  one  of  our  students,  Tim  Campbell,  in  collaboration  with  Dr.  Aiichiro  Nakano  in 
the  Department  of  Computer  Science.  In  this  scheme,  the  computational  workload  is  distributed 
equally  among  the  processors  and  the  interprocessor  communication  is  minimal.  In  this  scheme, 
the  atomic  coordinates  are  transformed  into  curvilinear  coordinates: 

^  =  X  -f-  X«Qexp[iQ-x)  .  (1) 

Q 

and  in  the  curvilinear  space  the  system  is  partitioned  uniformly  into  subsystems  which  are  mapped 
onto  the  processors  of  a  parallel  machine.  Using  the  simulated  annealing  approach,  the  parameters 
{a  }  are  varied  so  that  processors  have  equal  workloads  and  also  the  interprocessor 
communication  is  minimal.  (When  structural  rearrangement  occurs  during  a  simulation,  the 
variational  parameters  are  updated  dynamically  to  adjust  the  boundaries  of  subsystems.)  Figure  13 
shows  the  perfonnance  of  this  dynamic  load-balancing  algorithm  on  36  nodes  of  the  in-house 
Digital  Alpha  system  on  two  Giga  switches.  Evidently  the  algorithm  performs  very  well  even  when 
the  size  of  the  system  increases  from  25  to  125  million  atoms.  Recently  this  load  balancer  has  been 
further  improved  by  replacing  the  plane-wave  basis  by  orthonormal  wavelets. 
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Figure  13:  Performance  of  the  adaptive 
curvilinear-coordinate  load  balancer  on  the 
Digital  Alpha  system  with  two  Giga  switches. 
Each  node  has  3.5  million  atoms.  Note  the 
increase  in  the  execution  time  is  only  10%  even 
though  the  size  of  the  system  increases  from  25 
million  to  125  million  atoms 


§3  INTERACTIONS/TRANSITIONS 
Government  Laboratories 

We  have  a  number  of  collaborations  with  scientists  at  government  laboratories:  i)  We  are  working 
with  Drs.  Joseph  Lichtenhan  and  Stephen  Rodgers  at  the  Phillips  Laboratory;  and  ii)  with  Dr. 
Jeremy  Broughton  and  his  research  group  at  the  Naval  Research  Laboratory.  We  are  also 
collaborating  with  a  number  of  scientists  at  Argonne  National  Laboratory:  Dr.  C.  Loong  at  the 
Intense  Pulsed  Neutron  Source  Division  and  Drs.  T.  Disz,  W.  Gropp,  and  R.  Stevens  in  the 
Mathematics  and  Computer  Science  Division. 

Industries 

In  collaboration  with  Drs.  James  Patterson.  Joseph  Manke,  Daniel  Curtis,  Benjamin  Dembart. 
and  Thomas  Wicks  at  Boeing  Corporation,  one  of  our  students.  Timothy  Campbell,  developed  a 
parallel  implementation  of  electromagnetic  scattering  calculations  using  the  Fast  Multipole  Method 
(FMM). 


Universities 

Recently  we  put  together  an  outstanding,  multidisciplinary  team  consisting  of:  i)  D.  Brenner 
(North  Carolina  State);  ii)  E.  Kaxiras  (Harvard  Univ.):  hi)  J.  Broughton  (NRL);  iv)  A.  Madhukar 
(Univ.  of  Southern  California);  v)  L.  Greengard  (Courant);  and  vi)  T.  Disz,  W.  Gropp,  and  R. 
Stevens  (Argonne).  One  of  our  goals  is  to  develop  a  materials  simulation  software  repository 
which  will  include:  i)  Multiresolution  schemes  to  carry  out  large-scale  MD  simulations  based  on 
interatomic  potentials  derived  from  electronic-stmcture  calculations;  ii)  new  integration  algorithms 
that  will  allow  MD  simulations  to  reach  microsecond  time  scales;  hi)  extensions  in  the  areas  of 
dynamic  process  support,  remote  memory  operations,  and  various  latency-hiding  operations  that 
will  allow  substanti^  progress  in  computations  on  heterogeneous  pai'allel  platforms;  and  iv) 
interactive  and  immersive  visualization  tools  for  materials  simulations. 

§4  EDUCATIONAL  ACCOMPLISHMENTS:  TRAINING  OF  GRADUATE 
STUDENTS 

Recently  we  introduced  an  interdisciplinary  program  which  allows  graduate  students  to  obtain  a 
Ph.D.  in  physics  and  a  M.S.  from  the  Department  of  Computer  Science.  The  aim  of  this  program 
is  to  provide  students  with  broad-based  training  in  high  performance  computing  and 
communications  (HPCC)  and  the  physical  sciences.  In  connection  with  this  program,  we  have 
introduced  a  number  of  high  performance  computing  courses  in  the  Physics  and  Computer  Science 
Departments.  The  Department  of  Physics  has  two  graduate  courses  in  computational  physics 
which  are  cross-listed  with  computer-science  courses.  The  first  course  deals  with  classical  and 
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quantum  simulations  on  parallel  architectures.  The  second  course,  designed  for  advanced  graduate 
students,  covers  special  topics  such  as  multiscale  phenomena,  multigrid  methods,  wavelets,  etc. 
In  the  Department  of  Computer  Science  we  have  introduced  three  new  HPCC  courses.  Two  more 
courses  will  soon  be  added  to  the  computer-science  curriculum:  a)  Heterogeneous  Computing;  and 
b)  Scientific  Visualization.  The  courses  we  have  introduced  emphasize  parallel  computing  and 
algorithm  design  for  large-scale  scientific  applications.  Students  have  access  to  a  number  of 
parallel  machines  to  gain  hands-on  experience  and  to  perform  research  on  large-scale  computational 
projects. 


Concurrent  Computing  Laboratory  for  Materials  Simulations 


MP 

6,1^^  nocl^s 
SiMt)  arcbiJtecfurte 


155  Mbps 


High  Speed  Network/ATM  ^  SUN  S670 


Figure  14:  Computational  infrastructure  at  our  Concurrent  Computing  Laboratory  for  Materials  Simulations 
(CCLMS).  The  CCLMS  consists  of  two  parallel  computing  laboratories,  one  in  Physics  and  Astronomy  and  the 
other  in  Computer  Science.  With  $2.5  million  in  infrastructure  enhancement  grants  from  the  State  of  Louisiana, 
these  labs  have  been  equipped  with  the  following  parallel  machines:  Digital  Alpha  cluster  -40  Alpha 
nodes  (46  processors  consisting  of  eight  275-MHz  and  thirty-eight  165-MHz  processors)  connected  via  two  Giga 
switches;  Intel  iWarp  -  a  64-cell  systolic  architecture;  MasPar-  an  8192-node  SIMD  (Single  Instruction 
Multiple  Data)  machine;  Intel  iPSC/860  -a  distributed-memory  MIMD  (Multiple  Instructions  Multiple 
Data)  machine.  The  CCLMS  also  has  visualization  facilities  which  include  a  Silicon  Graphics  workstation,  an  8- 
processor  Silicon  Graphics  Power  Center,  and  stereo-graphic  equipment.  Efforts  are  underway  to  establish  a  virtual 
environment  (VE)  laboratory  which  will  feature  an  interactive  and  immersive  ImmersaPesk  for  visualization. _ 


Our  students  also  have  excellent  opportunities  to  broaden  their  research  experience  beyond  the 
traditional  university  based  environment.  They  are  involved  in  our  collaborative  efforts  with 
computational  and  experimental  scientists  at  government  laboratories,  industries,  and  other 
universities.  These  interactions  have  significantly  enhanced  the  research  capabilities  of  students. 
Through  these  contacts,  students  also  have  access  to  excellent  parallel  computing  and  visualization 
facilities  at  other  institutions. 
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